d <- readRDS(file='Rdata/input/data_2020-04-12.rds')
#
country <- as.character(d$Country)
# mean sd of ditrivution infection to death
mu_id = 18.8
si_std = .45*mu_id
load and clean mobility data.
## [1] "Algeria" "China"
## [3] "Czechia" "Dominican_Republic"
## [5] "Ecuador" "Iran"
## [7] "Peru" "South_Korea"
## [9] "United_Kingdom" "United_States_of_America"
## [1] "warning"
## [1] "Algeria" "China" "Dominican_Republic"
## [4] "Ecuador" "Iran" "Peru"
## [7] "South_Korea"
# correct epi
D <- d$D_active_transmission
mD <- as.matrix(D[,-1])
I <- d$I_active_transmission
cumD <- data.frame(D$dates,apply(D[,-1],2,cumsum))
cumI <- data.frame(I$dates,apply(I[,-1],2,cumsum))
N_geo <- ncol(D)-1
N_days <- nrow(D)
country <- as.character(d$Country)
load data (dots), get weekly average, interpolate between (blue line)
load Epi data and serial interval
## [1] 0.0292052
see paper
# delay death infection matrix
H <- matrix(0,nrow = N_days, ncol = N_days)
for (i in 1:N_days){
f <- max(c(1,i-delta_id$SItrunc))
H[i,f:i] <- rev(delta_id$dist)[((delta_id$SItrunc+1)-(i-f)):(delta_id$SItrunc+1)]
if (i>1) H[i,f:i] <- H[i,f:i]/sum( H[i,f:i])
}
# delay SI death to death matrix
W <- matrix(0,nrow = N_days, ncol = N_days)
for (i in 1:N_days){
f <- max(c(1,i-SI$SItrunc))
W[i,f:i] <- rev(SI$dist)[((SI$SItrunc+1)-(i-f)):(SI$SItrunc+1)]
# if (i>1) H[i,f:i] <- H[i,f:i]/sum( H[i,f:i])
}
# overall infectivity matrix
Ot <- W%*%mD
# correct infectivity when >0 case observed but no infectivity
for(i in 1:N_geo){
Ot[ which( (D[,i+1]>0) & (Ot[,i] ==0 ) ),i] <-NA
}
# matrix of mobility and delayed version
M <- as.matrix(mobility$mob_combined_smooth[,-1])
# M_D <- H %*% M
get epiestim Rt estimates with 7 day time window.
# epi
R0 <- 5
# for risk
b <- 0.5
theta0 <- c(rep(R0,N_geo),
rep(b,N_geo))
prior_theta <- matrix(c(rep(c(0,10),N_geo),
rep(c(0,1e3),N_geo)),
length(theta0),2, byrow=TRUE)
# parameter names
f0 <- function(x) paste0('R0_',x)
f1 <- function(x) paste0('beta_',x)
n_t<- c(sapply(country,f0), sapply(country,f1))
# sd dev for proposal
sigma <- rep(1e-1,length(theta0))
# useful functions
sapply(paste0('Rscript/MCMC/',(list.files('Rscript/MCMC/'))),FUN = source)
## Rscript/MCMC/adapt.R Rscript/MCMC/adapt_tuning.R
## value ? ?
## visible FALSE FALSE
## Rscript/MCMC/Like1.R Rscript/MCMC/MCMC_full.R
## value ? ?
## visible FALSE FALSE
## Rscript/MCMC/MCMC_iter.R
## value ?
## visible FALSE
#check
rep <- 2e4
# res <- MCMC_iter(iter = rep, theta0 = theta0, s = sigma)
res <- MCMC_full(iter = rep, theta0 = theta0, s = sigma, repli_adapt = 10, within_iter = rep/10)
# save.image(file = 'Rdata/check_mobility.RData')
# load(file = 'Rdata/check_mobility.RData')
Acc <- colSums(diff(res$theta)!=0)/(rep-1)
plot(res$logL[,1])
layout(matrix(1:4,2,2))
for (i in 1:length(theta0)){
plot(res$theta[,i],
main = paste0(n_t[i],' - ',round(Acc[i]*100)))#, ylim = prior_theta[i,])
}
Rt_daily <- array(NA,dim = c(N_days,N_geo,rep))
Rt_D2 <- Rt_daily
m_eff <- Rt_daily
for (j in 1:rep){
R_daily <- Rt_fun(theta = res$theta[j,], x = M )
Rt_daily[,,j] <- R_daily
Rt_D2[,,j] <- H %*% R_daily
B <- rep(1,nrow(M)) %*% t(res$theta[j,(2-1)*N_geo+ (1:N_geo)])
m_eff[,,j] <- 1 + 1/B * log( H %*% exp( - B * (1-M) ) )
}
n_d <- 1e2
Rt_assumed_mob <- array(NA,dim = c(n_d,N_geo,rep))
x <- matrix(seq(0,1,length.out = n_d),n_d,N_geo,byrow = FALSE)
for (j in 1:rep){
R_daily <- Rt_fun(theta = res$theta[j,], x = 1-x )
Rt_assumed_mob[,,j] <- R_daily
}
# format outputs
temp <- D
temp[,-1] <- NA
results_full_Rt_daily <- list(median_R = temp,
low_R = temp,
up_R = temp,
M = temp)
results_full_Rt_D2 <- list(median_R = temp,
low_R = temp,
up_R = temp,
M_D = temp)
results_meff <- list(median_meff = temp,
low_meff = temp,
up_meff = temp)
for (i in 1:N_geo){
temp <- apply(Rt_daily[,i,],1,quantile,c(.5,.025,.975),na.rm = TRUE)
results_full_Rt_daily$median_R[,i+1] <- temp[1,]
results_full_Rt_daily$low_R[,i+1] <- temp[2,]
results_full_Rt_daily$up_R[,i+1] <- temp[3,]
results_full_Rt_daily$M[,i+1] <- M[,i]
temp <- apply(Rt_D2[,i,],1,quantile,c(.5,.025,.975),na.rm = TRUE)
results_full_Rt_D2$median_R[,i+1] <- temp[1,]
results_full_Rt_D2$low_R[,i+1] <- temp[2,]
results_full_Rt_D2$up_R[,i+1] <- temp[3,]
# results_full_Rt_D2$M_D[,i+1] <- M_D[,i]
temp <- apply(m_eff[,i,],1,quantile,c(.5,.025,.975),na.rm = TRUE)
results_meff$median_meff[,i+1] <- temp[1,]
results_meff$low_meff[,i+1] <- temp[2,]
results_meff$up_meff[,i+1] <- temp[3,]
}
# format outputs
temp <- data.frame(matrix(NA,nrow = n_d, ncol = 1+N_geo))
temp[,1] <- x[,1]
names(temp) <- c('mobility', country)
results_full_Rt_assumed_mob <- list(median_R = temp,
low_R = temp,
up_R = temp)
for (i in 1:N_geo){
temp <- apply(Rt_assumed_mob[,i,],1,quantile,c(.5,.025,.975),na.rm = TRUE)
results_full_Rt_assumed_mob$median_R[,i+1] <- temp[1,]
results_full_Rt_assumed_mob$low_R[,i+1] <- temp[2,]
results_full_Rt_assumed_mob$up_R[,i+1] <- temp[3,]
}
n_pred <- 7
rep_sim <- 1e2
# choose future mobility
M_pred <- rbind(M,matrix(M[nrow(M),],n_pred,N_geo,byrow = TRUE))
# useful functions
sapply(paste0('Rscript/projections/',(list.files('Rscript/projections/'))),FUN = source)
## Rscript/projections/mobility_prediction.R
## value ?
## visible FALSE
res_projection <- mobility_prediction(n_pred = n_pred, rep_sim = rep_sim)
temp <- data.frame(dates = seq(1,n_pred) + D$dates[nrow(D)],
matrix(NA,n_pred,N_geo))
names(temp) <- c('dates',country)
summary_proj <- list(median = temp,
low = temp,
up = temp)
for (i in 1:N_geo){
temp <- apply(res_projection$D_pred[,i,],1,quantile,c(.5,.025,.975))
summary_proj$median[,i+1] <- temp[1,]
summary_proj$low[,i+1] <- temp[2,]
summary_proj$up[,i+1] <- temp[3,]
}
n_pred <- 60
rep_sim <- 1e2
# choose future mobility
M_pred <- rbind(M,matrix(0.05,n_pred,N_geo,byrow = TRUE))
res_projection_LT1 <- mobility_prediction(n_pred = n_pred, rep_sim = rep_sim)
temp <- data.frame(dates = seq(1,n_pred) + D$dates[nrow(D)],
matrix(NA,n_pred,N_geo))
names(temp) <- c('dates',country)
summary_proj_LT1 <- list(median = temp,
low = temp,
up = temp)
for (i in 1:N_geo){
temp <- apply(res_projection_LT1$D_pred[,i,],1,quantile,c(.5,.025,.975))
summary_proj_LT1$median[,i+1] <- temp[1,]
summary_proj_LT1$low[,i+1] <- temp[2,]
summary_proj_LT1$up[,i+1] <- temp[3,]
}
n_pred <- 60
rep_sim <- 1e2
# choose future mobility
new_M <- matrix(NA,n_pred,N_geo,byrow = TRUE)
for (i in 1:N_geo){
x <- 0:(n_pred-1)
new_M[,i] <- (1-M[nrow(M),i])/(n_pred-1) * x + M[nrow(M),i]
}
M_pred <- rbind(M,new_M)
res_projection_LT2 <- mobility_prediction(n_pred = n_pred, rep_sim = rep_sim)
temp <- data.frame(dates = seq(1,n_pred) + D$dates[nrow(D)],
matrix(NA,n_pred,N_geo))
names(temp) <- c('dates',country)
summary_proj_LT2 <- list(median = temp,
low = temp,
up = temp)
for (i in 1:N_geo){
temp <- apply(res_projection_LT2$D_pred[,i,],1,quantile,c(.5,.025,.975))
summary_proj_LT2$median[,i+1] <- temp[1,]
summary_proj_LT2$low[,i+1] <- temp[2,]
summary_proj_LT2$up[,i+1] <- temp[3,]
}
for (i in 1:N_geo){
niceplot(i)
}
temp <- summary_proj
layout(matrix(1:4,2,2))
for (i in 1:N_geo){
plot_proj(i,threshold = FALSE,log = FALSE)
}
temp <- summary_proj_LT1
layout(matrix(1:4,2,2))
for (i in 1:N_geo){
plot_proj(i,threshold = TRUE,log = FALSE)
}
temp <- summary_proj_LT2
layout(matrix(1:4,2,2))
for (i in 1:N_geo){
plot_proj(i,threshold = FALSE,log = TRUE)
}
f <- which(country %in% 'United_Kingdom')
layout(matrix(1:4,2,2))
temp <- summary_proj
plot_proj(f,threshold = FALSE,log = FALSE)
temp <- summary_proj_LT1
plot_proj(f,threshold = TRUE,log = FALSE)
temp <- summary_proj_LT2
plot_proj(f,threshold = TRUE,log = TRUE)
when contained??
contained <- data.frame(country = country,
median = rep(NA,N_geo),
low = rep(NA,N_geo),
up = rep(NA,N_geo))
for (i in 1:N_geo){
f <- c(tail( which(results_full_Rt_assumed_mob$median_R[i+1]>1) , 1 ),
tail( which(results_full_Rt_assumed_mob$median_R[i+1]>1),1) +1 )
contained$median[i] <- mean(results_full_Rt_assumed_mob$median_R[f,1])
f <- c(tail( which(results_full_Rt_assumed_mob$low_R[i+1]>1) , 1 ),
tail( which(results_full_Rt_assumed_mob$low_R[i+1]>1),1) +1 )
contained$low[i] <- mean(results_full_Rt_assumed_mob$low_R[f,1])
f <- c(tail( which(results_full_Rt_assumed_mob$up_R[i+1]>1) , 1 ),
tail( which(results_full_Rt_assumed_mob$up_R[i+1]>1),1) +1 )
contained$up[i] <- mean(results_full_Rt_assumed_mob$up_R[f,1])
}
write.csv(x = contained,file = 'Rdata/results_contained.csv')
save.image(file = 'Rdata/check_mobility.RData')